Flow correlated percolation during vascular network formation in tumors 
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A theoretical model based on the molecular interactions between a growing tumor and a dynamically evolv- 
ing blood vessel network describes the transformation of the regular vasculature in normal tissues into a highly 
inhomogeneous tumor specific capillary network. The emerging morphology, characterized by the compartmen- 
talization of the tumor into several regions differing in vessel density, diameter and necrosis, is in accordance 
with experimental data for human melanoma. Vessel collapse due to a combination of severely reduced blood 
flow and solid stress exerted by the tumor, leads to a correlated percolation process that is driven towards criti- 
cality by the mechanism of hydrodynamic vessel stabilization. 

PACS numbers: 87.18.-h, 87.10+e, 87.17.Aa, 61.43Hv 



Tumor vasculature, the network of blood vessels in and 
around a growing tumor, is in many respects different from 
the regular vasculature in normal tissues. Hypoxia, the lack 
of oxygen, that prevents a small tumor nucleus from further 
growth, induces the expression of various diffusible growth 
factors (GF) by the tumor cells that trigger a coordinated re- 
sponse of angiogenesis - the formation of irregular blood ves- 
sels (for a review see |1, 2]). The expected increase in mi- 
crovasular density (MVD) is usually observed in the periph- 
ery of the tumor, whereas the morphology of the vasculature 
in the central part of the tumor is characterized by a decreased 
MVD, dilated vessels and regions of necrotic tumor tissue 
fs^ |4i] . The resulting tumor specific capillary network is very 
heterogeneous, composed of dense and void regions, and has a 
fractal dimension different from normal arteriovenous or nor- 
mal capillary networks ||5|]. 

Although on the molecular level the main actors in the an- 
giogenic game are rapidly identified, the physical principles 
that determine the global morphology of the vascular network 
in tumor tissues are not known. Since for instance MVD is 
used as a diagnostic tool in cancer therapy |6] a quantitative 
understanding of the mechanism that leads to the compart- 
mentalization of the tumor vasculature into various regions 
differing substantially in vessel density appears mandatory. 
Moreover, scale-invariant aspects like fractal dimension, are 
used as hints towards the nature of the growth process un- 
derlying the formation of the tumor vasculature |7]. In this 
Letter we propose a theoretical model for the evolution of tu- 
mor vasculature that illuminates the physical principles lead- 
ing to its global morphology. The experimentally observed 
increase in MVD at the tumor perimeter and periphery and 
decrease in MVD and vessel dilation in the tumor center in 
human melanoma |4] appears also as the general scenario in 
the theoretical model that we discuss. Furthermore, we will 
argue that vessel collapses in the interior of the tumor lead to 
a percolation process which is driven towards criticality, the 
percolation threshold, via a mechanism of vessel stabilization 
by increased blood flow in the remaining vessels. 



Guided by a two-dimensional cellular automaton model 
that two of us developed recently |8] we consider the tumor- 
vessel system as a dynamically evolving network or graph in- 
teracting with a tumor growth process (inspired by the Eden 
model 1 9]). The interaction takes place via two concentration 
fields: the oxygen originating in the vessel network, and the 
growth factor originating in the tumor cells (TC). A hydro- 
dynamic flow is imprinted on the vessel network that emits 
oxygen. TC's proliferate/die when the local oxygen concen- 
tration is high/low. Vessels (edges) emerge when the local 
GF concentration is high enough, and they vanish (collapse) 
stochastically inside the tumor, when the hydrodynamic shear 
force acting on the vessel walls is too low. The biological and 
pathophysiological motivation for the details of the model def- 
inition to follow is discussed in 1 8]. 

To be specific, we describe the topology of the vessel net- 
work by a graph G = (V, and identify each edge e^E with 
a vessel, and each node v G V with a vessel junction, where 
more than two vessels merge. Here we restrict to capillary net- 
works and do not discriminate between arteries and veins. The 
network G is embedded in the three-dimensional Euclidean 
space and restricted to the cube ZdR^ of volume L ^ . This 
cube is discretized into = {l /af" unit cells, where a dimen- 
sionless vector r = xi + + zk with x, z = 0, 1 , 2, . . . ,L — 1 
denotes each unit cell. The microscopic length scale is chosen 
to be (2 = lO^m, the typical size of the endothelial cell (EC) 
and TC. 

For computational convenience we restrict the edges to 
run only parallel to the three coordinate axes and identify 
an edge with the string of unit cells of Z that it covers: Let 
r^(^) and r^(^) be the two end-points of an edge e and 
l{e) — \vt(e) — Vs(e) \ =n(e) — \ the length of the vessel, then 

^ = {r = r(r,(^),r,(^), Q | C = 0, 1,2, . . . ,£(^)} (1) 

with r(r„r,,Q = r, + ^(r, - r,)/^(^). Note that V = 
{r,(^)|^G^}U{r,(^)|^G^}. 

The tumor is represented by the set T of points that are 
occupied by tumor cells: T = {r|A TC exists at r}. The ves- 
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sel network G is the source of an oxygen concentration field 
O2 (r) and the tumor T is the source of a growth factor con- 
centration field GF{r): 

GF{r) = £/zs^^(|r-r'|), 
r'er 

02(r) = I/z«oxy(|r-r'|). 

r'eE 

Rgf and 7?oxy are the growth factor and oxygen diffusion radii, 
respectively, and for simplicity we choose a piecewise linear 
and normalized form for the contribution /z/^(r) of each tumor 
cell / vessel segment, /z/?(r) = (1 — r / R) / {nR^ / 3) for r < 7? 
and hR{r) = for r > 7?, satisfying Jq drhR{r)47lR^ = 1. 

Each edge e represents a tubular vessel of diameter d{e), 
through which a hydrodynamic blood flow of magnitude q{e), 
exerting a shear force f{e) upon the vessel walls, can pass. 
The flow is assumed to be an incompressible laminar station- 
ary flow for which q{e) and /(^) are calculated by Poiseuille's 
law: 

q{e)=d^{e)VP{e), and f{e)=d{e)VP{e), (2) 

where the pressure gradient in the vessel e is defined via 
yp{e) = \P{rt{e)) - P{rs{e))\ and the pressure P{r) in the 
nodes (vessel junctions) of the network is computed using 
Kirchhoff 's law. The boundary condition for the blood pres- 
sure P{r) is static and chosen in such a way that blood flow as 
well as shear force in the original network is homogeneous: 
P{r = {x,y,z)) = 1.5|x + j + z|/[3(L- 1)] at the boundary 
3Z = {r = {x^y^z)\x = or 3; = or z = 0}. 

Initial configuration: The original tissue is regularly vascu- 
larized with a homogeneous capillary network of given MVD 
that is fixed by inter- vessel distance 5: E = Eq = {e\rs{e) = 
b{n\,n2,m),Yt{e) = r^(^) + 5i,5j, or 5k, ^1,^2,^3 = 
0, 1, . . . - 1} with = [L/5J + 1. The number of nodes 
n{V) isN^, that of edges is 3(A^- 1)A^^ and^(^) =5+1. 
For each edge e, d{e) = l(10/im), q{e) = qo, and f{e) = fo 
with ^0 = /o = 0.5 /L from Eq. A tumor nucleus contain- 
ing Ntc tumor cells and grown using the Eden rule 1 9] starting 
with a seed at the system center = (L/2,L/2,L/2) defines 
the set of tumor cells T at time ^ = 0. ^uo(r) = for all r G T, 
representing the time spent in hypoxia. Starting with this ini- 
tial configuration the following computations are performed 
sequentially in each time step of duration At = Ih. 

TC proliferation: Proliferation of tumor cells is possible at 
tumor surface sites r G 5 = {r|r ^ T, r' G T, |r' — r| = 1} if the 
local O2 concentration is sufficiently large: T ^TU {r} with 
probability At/tjc if ^2(r) > <^oxy. 

Vessel growth: New vessels (edges in G) can be grown 
(added to G) in regions, where the local GF concentration is 
sufficiently large: /) Randomly choose ei e E and ri e ei. 
a) If there is r2 G ^2 for a certain ^2 7^ ^1 so that r2 — ri || i 
or j or k, 2 < |r2 -ri| < ^max, GF(r(ri,r2, Q) > Cgf for aU 
^ = 0, 1, . . . , |r2 — ri|, and r(ri,r2,Q ^ e for any e e E with 
^ = 1, . . . , |r2 - ri I - 1, then E ^ EU {^new} with ^new 
{r|r = r(ri,r2,Q,C = 0, 1,2, . . . , |r2 - ri |}. Repeat /) and //) 



3n{E)At /tEc times, since the number of potential sprouting 
events should be proportional to the total length of the vessels 
in the network. 

Vessel dilatation: Within the tumor no new vessels can oc- 
cur because of the lack of space, but vessels increase their 
diameter due to proliferation of EC's in the vessel walls if 
the local GF concentration is sufficiently large: For all e G 
E, d{e) d{e)^l^,^,e,t,^{GF{r)-Cgf)/{i{e) + l)/27l with 
probability At/tEc if ^(^) < <^max- Here 0step(-^) = 1 for x > 
and 0step (x) = for X < 0. 

Vessel regression and collapse: Weakly perfused vessels 
can collapse inside the tumor due to the solid stress ex- 
erted by the tumor: Compute P{r) for all r G V and then 
f{e) and q{e) according to Eq. (|2|). Vessels that are cut 
from the blood circulation (q{e) = 0) are instantaneously re- 
moved. For all other ^G£^weset£^^£^ — {^} with proba- 
bility A^/^coUapse if /(^)//o < ^Ic (i.e a weak shear force), and 
^({r|r G r, |r - r'| < 1 for r' G e}) > O.Si{e) (i.e. a large 
density of TCs around the vessel). 

TC death: TCs that are under-oxygenated longer than a 
time ^max will die: For all r G T, ^uo(r) ^uo(r) + 1 if ^2(1*) < 
Coxy. If ^uo(r) > Wx, T -{r} with probability 1/2. 

A schematic illustration of these procedures is shown in 
Fig. [1] We have simulated the model using various param- 
eter values, but here we restrict ourselves to the discussion of 
one typical parameter set, which is partly guided by data for 
human melanoma |4]. The original MVD is set by 5 = lOOjum 
and the oxygen diffusion radius Roxy = 100/im, yielding an 
average O2 concentration of O2 ~ 0.03. The oxygen thresh- 
old for TC proliferation and for cell death due to hypoxia is 
^^oxy = 0.01, i.e. clearly below O2. TC and EC proliferation 
time are ^tc = 10 h and tEC = 40 h, respectively, maximum TC 
survival time in hypoxia ^max = 20 h. The GF diffusion radius 
is Rgf = 200jLim, the GF threshold Cgf = 0.001. The maximum 
vessel diameter is Jmax = 35/im, maximum sprout migration 
distance ^max = lOO^m, critical shear force r|c = 0.5 and ves- 
sel collapse rate 1/Wiapse = 1/50 h. The initial tumor size 
A^TC = 27000 (i.e. an initial tumor diameter of ca. 0.6mm). 

An example for the time evolution of the tumor/vessel sys- 
tem in this model is shown in Fig. [l] Starting from a regu- 
lar vessel network the MVD in the peritumoral region is in- 
creased due to the supply of GFs from the tumor, as can best 
be seen in the snapshots of an equatorial cross section through 
the tumor center (last row of Fig.l). Once the tumor grows 
over this highly vascularized region, vessels start to collapse, 
by which the MVD in the interior of the tumor is continu- 
ously decreased until only a few thick vessels, surrounded by 
cuffs of TCs remain. Due to the reduced MVD in the tumor 
center regions become hypoxic and TCs will die leaving large 
necrotic regions. This compartmentalization of the tumor into 
different shells that can be discriminated by MVD, vessel di- 
ameter and necrosis is also observed in real tumors |4]. 

Fig.l shows that the original vascular network, consisting 
of capillaries of equal diameter arranged in a regular grid with 
a given MVD that guarantees homogeneous distribution of O2 
and a constant shear stress in all vessels, is dynamically trans- 
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FIG. 1: (a) Schematic illustration of the presented model: The upper 
panel indicates TC proliferation and death, the lower panel shows 
vessel formation, dilatation and regression/collpase (the modified 
part of the configuration is always indicated by broken lines), (b) 
The time evolution of the tumor- vessel system is demonstrated by 3 
snapshots at time t=0, 200, and 400. The upper panel shows only the 
tumor (note the necrotic regions inside), the middle panel only the 
vessel network (note the increased MVD at the tumor periphery, and 
the reduced MVD and dilated vessels in the tumor center), and the 
lower panel shows an equatorial cross section of the whole system 
in the xy plane at z = L/2 (lower). The parameters are chosen as 
mentioned in the text. The color code of the TCs represents the age 
scaled to [0, 1] and the color code of the vessel indicate the scaled 
blood flow, q{e)/qQ. 

formed into a compartmentalized network with irregularly 
arranged dilated vessels and characterized by an inhomoge- 
neous MVD and O2 distribution. This remodeling is strongly 
correlated with the blood flow pattern: When new vessels are 
generated, they share the blood flow with their parents, which 
causes all of them to have weaker shear forces, subject to po- 
tential vessel collapse. On the other hand, when such critical 
vessels are indeed removed, the blood flow has again to be 
redirected into the remaining vessels, leading to an increase 
of the shear forces in the involved vessels. At the same time, 
surviving vessels may increase their diameters, also resulting 
in higher shear forces. The ratio ^ec /^collapse basically controls 
this flow-correlated remodeling process consisting of genera- 
tion / dilatation and collapse and also affects the tumor vol- 
ume or necrosis: Necrosis dominates the tumor tissue in lack 




FIG. 2: Plots of the tumor density pxc^ MVD, vessel diameter d, 
pressure gradient VP, and shear force / as functions of the distance 
to the center P = |r — rc | and time t. The lower left plot shows the de- 
pendence of the pressure gradient on the azimuthal angle 9 between 
r — Tc and the diagonal — Fc. 

of oxygen supply for a larger value of the ratio while the total 
necrotic region shrinks for a smaller value. 

The dynamical evolution described thus far can be an- 
alyzed quantitatively by studying the following quanti- 
ties: The radial tumor density pTc(^) = ^{Tr) / n{ZR) , 
vessel density MVD(R) = n{ER)/{4nR^), vessel diame- 
ter d{R) = l,eeER^{e)/lleeER^ pressure gradient VP(7?) = 
I.eeER^P{e)/i:eeER^ and VP(e) = LeeE.^Pie) /HeeE,^ and 
shear force f{R) = Y,eeERf{^)/lleeER- Her we have used 
ZR = {r\R<\v-v,\<R^AR},ER = {eeE\R<\{Vs{e)^ 
rt{e))/2 - rd < R^ AR}, Eq = {e e E\e < cos'' f^^^ < 

0+Ae},andr/? = rnZ/?. 

As seen in Fig. |2l the peak of MVD{R) is in accordance 
with the tumor boundary 7?tc(0- The tumor grows approxi- 
mately spherically with Rjcit) — ^tc(O) — where the 
factor 2 is typical for the Eden growth. In the tumor cen- 
ter, MVD{R) and pTc(^) are both very low. Their shape 
similarity is due to continuous interactions between the tu- 
mor and the vessel network. Vessels that have long been ex- 
posed to GF produced by TCs have large diameters: d{R) in- 
creases linearly from 1 at 7? :^ 7?tc + ^gf to Jmax at the tumor 
center. Such a characteristic vessel morphology is also in a 
quantitative agreement with experimental data from the hu- 
man melanoma |4]. 

The blood pressure gradient in the tumor center is up to 
50% lower than in normal vessels which is, from hydrody- 
namic considerations, an immediate consequence of the in- 
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FIG. 3: Plots of A/g of Eq. J3| in the vessel networks at t = 400 for 
different values of ^collapse ^nd fmax with other parameters taking the 
typical values. The tumor and the peritumoral region extend up to 
R< 145, for which is computed. In all cases, the data are fitted as 
A/e ^ £-2.52(5) has a slope -2.52. The inset shows A^e 

measured in different shells of the same thickness 20 for ^collapse = 20 
and ^max = 20. The slopes of the upper dashed line and lower solid 
line are —2.24 and —1.68, respectively. 

crease MVD in the peritumoral region. Moreover, the pres- 
sure gradient is lowest in the direction orthogonal to the global 
flow (0 = 7r/2), since for a given network the flow tends to 
use shortest paths from source to sink. Finally, the shear force 
acting on each vessel wall depends on the vessel diameter and 
the pressure gradient. The sharp decrease of VP in the peri- 
tumoral region is inherited to f{R) while the increased vessel 
diameter in the tumor center gives rise to large shear forces in 
large vessels. 

The geometrical features of the emerging tumor vascula- 
ture in our model are obviously very different from the origi- 
nal, regular capillary network: It consist of a combination of 
dense and void regions that might possess fractal properties. 
We used the box-counting method to determine the fractal di- 
mension Df SiS 

Z)/ = -linilnA^e/lne, (3) 

where is the number of boxes of volume £^ necessary to 
cover the tumor vessel network, that is defined to lie within 
the outer limit of the peritumoral region. The plot of versus 
8 shown in Fig. |3l yields Z)/ = 2.52(5),in agreement with the 
value for the critical percolation cluster in the random bond- 
percolation process in three dimension 1 10]. We checked that 
we obtain the same value for a wide range of parameter values. 

From this observation we can conclude that the basic mech- 
anism responsible for the fractal properties of the tumor vas- 
culature in our model is the stochastic removal of vessels via 
vessel collapse and regression. In conventional percolation a 
critical cluster only emerges if the edge concentration is fixed 
to be exactly at the percolation threshold. In our model this 



fine tuning is not necessary: the dynamically evolving net- 
work drives itself into this critical state since the removal of 
vessels is correlated with the blood flow. The collapse of crit- 
ical vessels stabilizes the remaining ones due to an increase in 
blood flow, as shown in our quantitative analysis. 

It has been suggested |5] that the origin of the fractal ar- 
chitecture of tumor vasculature might be due to an underly- 
ing invasion percolation process 1 1 1 ] due to inhomogeneities 
in the growth supporting matrix. In view of the theoretical 
model we have presented, which does not involve any such 
matrix-inhomogeneities, we propose that it is rather the flow 
correlated percolation process due to collapsing vessels inside 
the tumor that determines the fractal properties of the tumor 
vasculature. A commonly accepted view, at least for a large 
class of tumors like melanoma, also shared by our theoreti- 
cal approach, is that neo-vascularization mainly occurs at the 
tumor perimeter and a drastic reduction of vessel density oc- 
curs in the interior of the tumor. In such a scenario it appears 
unlikely that the fractal properties attained during growth in 
the periphery, independent of having characteristics of inva- 
sion percolation or not, survive the random dilution process in 
the tumor center. Thus a theoretical model that does not take 
into account vessel collapse can possibly explain the observed 
fractal dimension of the vessel network in tumors where the 
central MVD is not drastically reduced but fails to do so in 
networks where it is. 

To conclude we have introduced a theoretical model for a 
dynamically evolving, three-dimensional vessel network in- 
teracting with a growing tumor, which is guided by exper- 
imental data for human melanoma. The emerging network 
morphology agrees well with those data and we find that the 
network is remodeled from a regular into a fractal structure 
with characteristics of random percolation. This suggests also 
for a large class of real solid tumor with decreased central 
MVD that the basic mechanism leading to the fractal features 
of the tumor vasculature is the mechanism of stochastic vessel 
collapse inside the tumor. 
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